# This file estimates effects of having ones first two children be of same sex for fathers

# install.packages("dplyr")
# install.packages("date")
# install.packages("AER")

library(dplyr)
library(date)
library(AER)

# Load data
load("fathers_firsttwo.rdata")

# Filter of if age of oldest is unkwown (select parents )
fathers <-
  fathers %>%
  filter(!is.na(age_oldest)) 

# split data file to each election 
# 2009 election
# Subset and create vote variable
# filter of parents with less than two children at the election and a firstborn older than 26

fathers09 <- 
  fathers %>%
  filter(! is.na(stemte_2009) ) %>%
  filter(age_oldest >= 5113 & no_children_09 > 1) %>%
  mutate(vote = stemte_2009 == "Stemte",
         no_children = no_children_09)

#create age dummies for fathers 
#dichotomize twin first
fathers09 <-
  within(fathers09, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })


fathers09 <-
  fathers09 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))


## Repeat for 13
fathers13 <- 
  fathers %>%
  filter(! is.na(stemt) ) %>%
  filter(age_oldest >= 6574 & no_children_13 > 1) %>%
  mutate(vote = stemt == "Stemte",
         no_children = no_children_13)

fathers13 <-
  within(fathers13, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

fathers13 <-
  fathers13 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

## Repeat for 14
fathers14 <- 
  fathers %>%
  filter(! is.na(ep_stemt) ) %>%
  filter(age_oldest >= 6574 & no_children_14 > 1) %>%
  mutate(vote = ep_stemt,
         no_children = no_children_14)

fathers14 <-
  within(fathers14, {
    age        <- ntile(date, 25)
    twin_first <- twin_first > 0
  })

fathers14 <-
  fathers14 %>%
  group_by(age) %>%
  mutate(child_age = ntile(age_oldest, 4))

# Run models 
# Simple models. 
mod09_r <- lm(vote ~ first_two_same, data = fathers09)
mod13_r <- lm(vote ~ first_two_same, data = fathers13)
mod14_r <- lm(vote ~ first_two_same, data = fathers14)

# Models with fixed effect for fathers and children's age
mod09 <- lm(vote ~ first_two_same + as.factor(age) * as.factor(child_age), data = fathers09)
mod13 <- lm(vote ~ first_two_same + as.factor(age) * as.factor(child_age), data = fathers13)
mod14 <- lm(vote ~ first_two_same + as.factor(age) * as.factor(child_age), data = fathers14)

## Instrument number of children with first child being twin
mod09_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     first_two_same + as.factor(age) * as.factor(child_age), data = fathers09)
mod13_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     first_two_same + as.factor(age) * as.factor(child_age), data = fathers13)
mod14_iv  <- ivreg(vote ~ no_children + as.factor(age) * as.factor(child_age) | 
                     first_two_same + as.factor(age) * as.factor(child_age), data = fathers14)

## Run first-stage models 
mod09_fs <- lm(no_children ~ first_two_same + as.factor(age) * as.factor(child_age), data = fathers09)
mod13_fs <- lm(no_children ~ first_two_same + as.factor(age) * as.factor(child_age), data = fathers13)
mod14_fs <- lm(no_children ~ first_two_same + as.factor(age) * as.factor(child_age), data = fathers14)


# Compile results 

table_dad_09 <-
  rbind(summary(mod09_r)$coefficients,
        summary(mod09)$coefficients[2,],
        summary(mod09_iv)$coefficients[2,])

table_dad_13 <-
  rbind(summary(mod13_r)$coefficients,
        summary(mod13)$coefficients[2,],
        summary(mod13_iv)$coefficients[2,])

table_dad_14 <-
  rbind(summary(mod14_r)$coefficients,
        summary(mod14)$coefficients[2,],
        summary(mod14_iv)$coefficients[2,])

tab_first_stage <- 
  rbind(summary(mod09_fs)$coefficients[1:2,],
        summary(mod13_fs)$coefficients[1:2,],
        summary(mod14_fs)$coefficients[1:2,])


outtab_dad <- 
  rbind(cbind(table_dad_09, with(fathers09, table(first_two_same))),
        cbind(table_dad_13, with(fathers13, table(first_two_same))),
        cbind(table_dad_14, with(fathers14, table(first_two_same))),
        cbind(tab_first_stage, with(fathers14, table(first_two_same))))

write.table(outtab_dad, file = "father_firsttwo_res.txt")

sum_iv09 <- summary(mod09_iv, diagnostics = TRUE)
sum_iv13 <- summary(mod13_iv, diagnostics = TRUE)
sum_iv14 <- summary(mod14_iv, diagnostics = TRUE)

f_diag_fathers <- 
  rbind(sum_iv09$diagnostics[1,],
        sum_iv13$diagnostics[1,],
        sum_iv14$diagnostics[1,])

write.table(f_diag_fathers, file = "f_stat_fathers_firsttwo.txt")
